rm(list = ls())

## Libraries

suppressPackageStartupMessages(library(tidyverse, warn.conflicts = FALSE))
library(readr, warn.conflicts = FALSE)
library(psych, warn.conflicts = FALSE)
library(survey, warn.conflicts = FALSE)
library(broom, warn.conflicts = FALSE)
library(stargazer, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)

set.seed(221186)

## Load data

survey <- read_csv("../source_data/PACER+Ideology+Experiment_October+21%2C+2021_03.23.csv") %>% 
  slice(3:nrow(.)) %>%
  mutate(StartDate = as.POSIXct(StartDate),
         EndDate = as.POSIXct(EndDate),
         TimeElapsed = as.numeric(EndDate - StartDate)/60,
         attn_count = as.numeric(attn_check == "Other") + as.numeric(isspRoGGrid_pre_3 == "Against") + as.numeric(ds_5 == "Not disgusting at all"),
         treat = case_when(is.na(saw_pce_pre) & saw_control == "TRUE" ~ "Control",
                           saw_pce_pre == "TRUE" & saw_crip == "TRUE" ~ "PCE/CRIP",
                           saw_pce_pre == "TRUE" & saw_scf == "TRUE" ~ "PCE/SCF",
                           saw_pce_pre == "TRUE" & saw_sif == "TRUE" ~ "PCE/SIF",
                           saw_pce_pre == "TRUE" & saw_urip == "TRUE" ~ "PCE/URIP",
                           is.na(saw_pce_pre) & saw_crip == "TRUE" ~ "CRIP",
                           is.na(saw_pce_pre) & saw_scf == "TRUE" ~ "SCF",
                           is.na(saw_pce_pre) & saw_sif == "TRUE" ~ "SIF",
                           is.na(saw_pce_pre) & saw_urip == "TRUE" ~ "URIP"
         ),
         treat = ifelse(is.na(treat),"PCE",treat),
         treat2 = case_when(group %in% c(1, 11) ~ "Control",
                            group %in% c(2, 12) ~ "PCE",
                            group %in% c(3) ~ "CRIP",
                            group %in% c(4) ~ "PCE/CRIP",
                            group %in% c(5) ~ "URIP",
                            group %in% c(6) ~ "PCE/URIP",
                            group %in% c(7) ~ "SCP",
                            group %in% c(8) ~ "PCE/SCP",
                            group %in% c(9) ~ "SIP",
                            group %in% c(10) ~ "PCE/SIP"),
         ideo_treat_binary = group %in% c(3,4,5,6,7,8,9,10),
         pce_treat_binary = group %in% c(2,12,4,6,8,10),
         crip_binary = group %in% c(3,4),
         urip_binary = group %in% c(5,6),
         scp_binary = group %in% c(7,8),
         sip_binary = group %in% c(9,10),
         group = as.numeric(group),
         
         # Covariates
         
         past_vote_ge = factor(past_vote_ge, levels = c("Conservative", "Labour", "Liberal Democrats", "Green", "Brexit Party", "Plaid Cymru", "Scottish National Party (SNP)", "Other", "Did not vote", "Was not eligible to vote")),
         age = 2021 - as.numeric(`Year of Birth`),
         
         # Recode NAs which are Don't know to Don't know
         isspRoGGrid_pre_3 = ifelse(is.na(isspRoGGrid_pre_3), "Don't know",isspRoGGrid_pre_3),
         ds_5 = ifelse(is.na(ds_5), "Don't know",ds_5),
         
         # Recode outcome variables
         
         ## Set up raw variables
         
         redistSelf_pre_1_raw = redistSelf_pre_1,
         redistSelf_post_1_raw = redistSelf_post_1,
         taxSpendSelf_pre_1_raw = taxSpendSelf_pre_1,
         taxSpendSelf_post_1_raw = taxSpendSelf_post_1,
         perceptionsPoor_pre_1_raw = perceptionsPoor_pre_1,
         perceptionsPoor_post_1_raw = perceptionsPoor_post_1,
         perceptionsPoor_pre_2_raw = perceptionsPoor_pre_2,
         perceptionsPoor_post_2_raw = perceptionsPoor_post_2,
         perceptionsPoor_pre_3_raw = perceptionsPoor_pre_3,
         perceptionsPoor_post_3_raw = perceptionsPoor_post_3,
         bhpsRoGGrid_pre_1_raw = bhpsRoGGrid_pre_1,
         bhpsRoGGrid_pre_2_raw = bhpsRoGGrid_pre_2,
         bhpsRoGGrid_pre_3_raw = bhpsRoGGrid_pre_3,
         bhpsRoGGrid_post_1_raw = bhpsRoGGrid_post_1,
         bhpsRoGGrid_post_2_raw = bhpsRoGGrid_post_2,
         bhpsRoGGrid_post_3_raw = bhpsRoGGrid_post_3,
         isspRoGGrid_pre_1_raw = isspRoGGrid_pre_1,
         isspRoGGrid_pre_2_raw = isspRoGGrid_pre_2,
         isspRoGGrid_post_1_raw = isspRoGGrid_post_1,
         isspRoGGrid_post_2_raw = isspRoGGrid_post_2,
         
         ## redistSelf
         redistSelf_pre_1 = abs(as.numeric(redistSelf_pre_1) - 10),
         redistSelf_post_1 = abs(as.numeric(redistSelf_post_1) - 10),
         redistSelf_diff = redistSelf_post_1 - redistSelf_pre_1,
         
         ## taxSpendSelf
         taxSpendSelf_pre_1 = as.numeric(taxSpendSelf_pre_1),
         taxSpendSelf_post_1 = as.numeric(taxSpendSelf_post_1),
         taxSpendSelf_diff = taxSpendSelf_post_1 - taxSpendSelf_pre_1,
         
         ## perceptionsPoor_1
         perceptionsPoor_pre_1 = as.numeric(factor(perceptionsPoor_pre_1, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         perceptionsPoor_post_1 = as.numeric(factor(perceptionsPoor_post_1, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         perceptionsPoor_diff_1 = perceptionsPoor_post_1 - perceptionsPoor_pre_1,
         
         ## perceptionsPoor_2
         perceptionsPoor_pre_2 = as.numeric(factor(perceptionsPoor_pre_2, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         perceptionsPoor_post_2 = as.numeric(factor(perceptionsPoor_post_2, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         perceptionsPoor_diff_2 = perceptionsPoor_post_2 - perceptionsPoor_pre_2,
         
         ## perceptionsPoor_3
         perceptionsPoor_pre_3 = as.numeric(factor(perceptionsPoor_pre_3, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         perceptionsPoor_post_3 = as.numeric(factor(perceptionsPoor_post_3, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         perceptionsPoor_diff_3 = perceptionsPoor_post_3 - perceptionsPoor_pre_3,
         
         ## bhpsRoGGrid_1
         bhpsRoGGrid_pre_1 = as.numeric(factor(bhpsRoGGrid_pre_1, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         bhpsRoGGrid_post_1 = as.numeric(factor(bhpsRoGGrid_post_1, levels = rev(c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree")))),
         bhpsRoGGrid_diff_1 = bhpsRoGGrid_post_1 - bhpsRoGGrid_pre_1,
         
         ## bhpsRoGGrid_2
         bhpsRoGGrid_pre_2 = as.numeric(factor(bhpsRoGGrid_pre_2, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         bhpsRoGGrid_post_2 = as.numeric(factor(bhpsRoGGrid_post_2, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         bhpsRoGGrid_diff_2 = bhpsRoGGrid_post_2 - bhpsRoGGrid_pre_2,
         
         ## bhpsRoGGrid_3
         bhpsRoGGrid_pre_3 = as.numeric(factor(bhpsRoGGrid_pre_3, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         bhpsRoGGrid_post_3 = as.numeric(factor(bhpsRoGGrid_post_3, levels = c("Strongly disagree", "Disagree", "Neither agree nor disagree", "Agree", "Strongly agree"))),
         bhpsRoGGrid_diff_3 = bhpsRoGGrid_post_3 - bhpsRoGGrid_pre_3,
         
         ## jobsForAll
         jobsForAll_pre = bhpsRoGGrid_pre_3,
         jobsForAll_post = bhpsRoGGrid_post_3,
         jobsForAll_diff = bhpsRoGGrid_diff_3,
         
         ## isspRoGGrid_1
         isspRoGGrid_pre_1 = as.numeric(factor(isspRoGGrid_pre_1, levels = c("Strongly against", "Against", "Neither in favour of nor against", "In favour of", "Strongly in favour of"))),
         isspRoGGrid_post_1 = as.numeric(factor(isspRoGGrid_post_1, levels = c("Strongly against", "Against", "Neither in favour of nor against", "In favour of", "Strongly in favour of"))),
         isspRoGGrid_diff_1 = isspRoGGrid_post_1 - isspRoGGrid_pre_1,
         
         ## isspRoGGrid_2
         isspRoGGrid_pre_2 = as.numeric(factor(isspRoGGrid_pre_2, levels = c("Strongly against", "Against", "Neither in favour of nor against", "In favour of", "Strongly in favour of"))),
         isspRoGGrid_post_2 = as.numeric(factor(isspRoGGrid_post_2, levels = c("Strongly against", "Against", "Neither in favour of nor against", "In favour of", "Strongly in favour of"))),
         isspRoGGrid_diff_2 = isspRoGGrid_post_2 - isspRoGGrid_pre_2,
         
         ## Experienced COVID loss/support
         experienced_loss  = !is.na(`PACER:eKnowLoss_1`) | !is.na(`PACER:eKnowLoss_2`) | 
           !is.na(`PACER:eKnowLoss_3`) | !is.na(`PACER:eKnowLoss_4`),
         experienced_receipt = !is.na(`PACER:eKnowReceived_1`) | !is.na(`PACER:eKnowReceived_2`) | 
           !is.na(`PACER:eKnowReceived_3`) | !is.na(`PACER:eKnowReceived_4`)
         
         ) %>%
  filter(StartDate >= as.POSIXct("2021-08-24 01:12:58 BST"),
         DistributionChannel != "test")

## PCA extraction

redistIndex_princomp <- survey %>%
  select(redistSelf_pre_1, perceptionsPoor_pre_1, 
         perceptionsPoor_pre_2, perceptionsPoor_pre_3) %>%
  na.omit() %>%
  princomp(cor= TRUE)

RoGIndex_princomp <- survey %>%
  select(jobsForAll_pre, bhpsRoGGrid_pre_1,
         bhpsRoGGrid_pre_2, isspRoGGrid_pre_1,
         isspRoGGrid_pre_2) %>%
  na.omit() %>%
  princomp(cor = TRUE)

# output loadings tables

redistLoadings <- as.matrix.data.frame(redistIndex_princomp$loadings)

rownames(redistLoadings) <- dimnames(redistIndex_princomp$loadings)[[1]]
rownames(redistLoadings) <- c("redistSelf", "noFaultOfTheirOwn", "tooManyHandouts", "dontDeserveHelp")
colnames(redistLoadings) <- paste0("PCA",1:4)

sink("../generated_tex/redistribution_pca_loadings.tex")
xtable::print.xtable(xtable::xtable(round(redistLoadings, 3)), floating = FALSE)
sink()


rogLoadings <- as.matrix.data.frame(RoGIndex_princomp$loadings)

rownames(rogLoadings) <- dimnames(RoGIndex_princomp$loadings)[[1]]
rownames(rogLoadings) <- c("jobsForAll", "privateEnterprise", "StateOwnership", "govFinanceNewJobs", "supportNewProducts")
colnames(rogLoadings) <- paste0("PCA",1:5)

sink("../generated_tex/rog_pca_loadings.tex")
xtable::print.xtable(xtable::xtable(round(rogLoadings, 3)), floating = FALSE)
sink()


# construct PCA measures

survey <- survey %>%
  mutate(redistIndex_pre = predict(redistIndex_princomp,
                                   newdata = data.frame(redistSelf_pre_1 = redistSelf_pre_1,
                                                        perceptionsPoor_pre_1 = perceptionsPoor_pre_1,
                                                        perceptionsPoor_pre_2 = perceptionsPoor_pre_2,
                                                        perceptionsPoor_pre_3 = perceptionsPoor_pre_3))[,1],
         redistIndex_post = predict(redistIndex_princomp,
                                   newdata = data.frame(redistSelf_pre_1 = redistSelf_post_1,
                                                        perceptionsPoor_pre_1 = perceptionsPoor_post_1,
                                                        perceptionsPoor_pre_2 = perceptionsPoor_post_2,
                                                        perceptionsPoor_pre_3 = perceptionsPoor_post_3))[,1],
         redistIndex_diff = redistIndex_post - redistIndex_pre,
         
         RoGIndex_pre = predict(RoGIndex_princomp,
                                   newdata = data.frame(jobsForAll_pre = jobsForAll_pre,
                                                        bhpsRoGGrid_pre_1 = bhpsRoGGrid_pre_1,
                                                        bhpsRoGGrid_pre_2 = bhpsRoGGrid_pre_2,
                                                        isspRoGGrid_pre_1 = isspRoGGrid_pre_1,
                                                        isspRoGGrid_pre_2 = isspRoGGrid_pre_2))[,1],
         RoGIndex_post = predict(RoGIndex_princomp,
                                    newdata = data.frame(jobsForAll_pre = jobsForAll_post,
                                                         bhpsRoGGrid_pre_1 = bhpsRoGGrid_post_1,
                                                         bhpsRoGGrid_pre_2 = bhpsRoGGrid_post_2,
                                                         isspRoGGrid_pre_1 = isspRoGGrid_post_1,
                                                         isspRoGGrid_pre_2 = isspRoGGrid_post_2))[,1],
         RoGIndex_diff = RoGIndex_post - RoGIndex_pre)

## Treatment effects

dv_vars_lag <- c("redistSelf_post_1", "taxSpendSelf_post_1", "jobsForAll_post", "RoGIndex_post", "redistIndex_post")

iv_vars_lag <- c("redistSelf_pre_1", "taxSpendSelf_pre_1", "jobsForAll_pre", "RoGIndex_pre", "redistIndex_pre")


models_out <- lapply(1:length(dv_vars_lag), 
       function(i) {
         
         treatment_model <- survey %>% 
           filter(term_status == "complete") %>%
           mutate(dv = scale(get(dv_vars_lag[i])),
                  iv = scale(get(iv_vars_lag[i]))) %>%
           lm(dv ~ crip_binary + urip_binary + scp_binary + sip_binary + iv, data = .)
         
         null_model <- survey %>% 
           filter(term_status == "complete") %>%
           mutate(dv = scale(get(dv_vars_lag[i])),
                  iv = scale(get(iv_vars_lag[i]))) %>%
           lm(dv ~ iv, data = .)
         
         f_test <- anova(null_model, treatment_model)
         
         treatment_model_coefs <- tidy(treatment_model, conf.int = TRUE)
         treatment_model_coefs$p_value_f_test <- as.numeric(na.omit(f_test$`Pr(>F)`))
         treatment_model_coefs$model <- dv_vars_lag[i]
         treatment_model_coefs  
         
         list(treatment_model = treatment_model,
              null_model = null_model, 
              f_test = f_test,
              treatment_model_coefs = treatment_model_coefs)
         
         }) 


models_out_coefs <- lapply(models_out, function(x){
  x[["treatment_model_coefs"]]
  }
) %>% 
  plyr::rbind.fill() %>%
  as_tibble()


# Calculate bh adjustment

f_test_p_vals <- models_out_coefs %>% filter(term == "(Intercept)") %>% .$p_value_f_test

f_test_p_vals_adjusted <- p.adjust(f_test_p_vals, method = "BH")

# Coefficient table

sink("../generated_tex/experiment_treatment_effects.tex")
stargazer(lapply(models_out, function(x) x$treatment_model),
          covariate.labels = c("Common Risk and Insurance",
                               "Unequal Risk and Insurance",
                               "State Capacity",
                               "State Incapacity",
                               "Pre-treatment DV"),
          add.lines = list(c("F-test adj. p-value", round(f_test_p_vals_adjusted,3))),
          keep.stat = c("N", "rsq"),
          no.space = TRUE,
          column.labels = gsub("_post_1|_post","",dv_vars_lag),
          dep.var.labels = "",
          omit.table.layout = "n",
          report = c("vcs"), 
          model.numbers = FALSE,
          font.size = "footnotesize",
          float = FALSE)

sink()

# Plot the coefficients

x <- models_out_coefs %>% 
  filter(term %in% c("crip_binaryTRUE", "urip_binaryTRUE", "scp_binaryTRUE", "sip_binaryTRUE")) %>%
  mutate(term2 = factor(case_when(term == "crip_binaryTRUE" ~ "Common Risk and Insurance",
                           term == "urip_binaryTRUE" ~ "Unequal Risk and Insurance",
                           term == "scp_binaryTRUE" ~ "State Capacity",
                           term == "sip_binaryTRUE" ~ "State Incapacity"),
                        levels = rev(c("Common Risk and Insurance", "Unequal Risk and Insurance", "State Capacity", "State Incapacity"))),
         model2 = factor(case_when(model == "jobsForAll_post" ~ "jobsForAll",
                            model == "redistIndex_post" ~ "redistIndex",
                            model == "redistSelf_post_1" ~ "redistSelf",
                            model == "RoGIndex_post" ~ "RoGIndex",
                            model == "taxSpendSelf_post_1" ~ "taxSpendSelf"),
                         levels = c("jobsForAll", "redistSelf", "taxSpendSelf", "redistIndex", "RoGIndex"))) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, group = model2, y = term2)) + 
    geom_pointrange() + 
    facet_wrap(~model2, ncol = 5) +
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  xlab("Estimate") + 
  ylab("")

ggsave("../generated_images/treatment_effects.png", plot = x, width = 10, height = 4)


## Treatment x PCE effects

models_out_int <- lapply(1:length(dv_vars_lag), 
                     function(i) {
                       
                       treatment_model <- survey %>% 
                         filter(term_status == "complete") %>%
                         mutate(dv = scale(get(dv_vars_lag[i])),
                                iv = scale(get(iv_vars_lag[i]))) %>%
                         lm(dv ~ ideo_treat_binary * pce_treat_binary + iv, data = .)
                       
                       
                       sims <- data.frame(mvtnorm::rmvnorm(2000, coef(treatment_model), vcov(treatment_model)))
      
                       PCConsistent <- (sims$ideo_treat_binaryTRUE + sims$ideo_treat_binaryTRUE.pce_treat_binaryTRUE) > 0 &
                         (sims$ideo_treat_binaryTRUE.pce_treat_binaryTRUE > 0)
                       
                       PCConsistent_pct <- mean(PCConsistent)
                       
                       p_pc_consistent <- 1 - PCConsistent_pct
                       
                       treatment_model_coefs <- tidy(treatment_model, conf.int = TRUE)
                       treatment_model_coefs$p_pc_consistent <- p_pc_consistent
                       treatment_model_coefs$model <- dv_vars_lag[i]
                       treatment_model_coefs  
                       
                       names(p_pc_consistent) <- dv_vars_lag[i]
                       
                       list(treatment_model = treatment_model,
                            p_pc_consistent = p_pc_consistent, 
                            treatment_model_coefs = treatment_model_coefs)
                       
                     }) 

# Calculate bh adjustment

p_vals_adjusted <- p.adjust(unlist(lapply(models_out_int, function(x) x$p_pc_consistent)), method = "BH")


# Coefficient table

sink("../generated_tex/experiment_treatment_interaction_effects.tex")
stargazer(lapply(models_out_int, function(x) x$treatment_model),
          covariate.labels = c("Ideological Treatment",
                               "Personal Experience Treatment",
                               "Pre-treatment DV",
                               "Interaction"),
          add.lines = list(c("Adj. p-value", round(p_vals_adjusted,3))),
          keep.stat = c("N", "rsq"),
          no.space = TRUE,
          column.labels = gsub("_post_1|_post","",dv_vars_lag),
          dep.var.labels = "",
          omit.table.layout = "n",
          report = c("vcs"), 
          model.numbers = FALSE,
          font.size = "footnotesize",
          float = FALSE)

sink()


x <- lapply(models_out_int, function(x) x$treatment_model_coefs ) %>%
  plyr::rbind.fill() %>% as_tibble() %>%
  filter(term == "ideo_treat_binaryTRUE:pce_treat_binaryTRUE") %>%
  mutate(model2 = gsub("_post_1|_post","",model),
         model2 = factor(model2, levels = model2[order(estimate)])) %>%
  ggplot(aes(x = estimate, y = model2, xmin = conf.low, xmax = conf.high)) + 
  geom_pointrange() + 
  theme_bw() + 
  xlab("Estimated ideology * personal experience effect") + 
  ylab("") + 
  geom_vline(xintercept = 0, linetype = 2)

ggsave("../generated_images/treatment_interaction_effects.png", plot = x, width = 6, height = 4)


## Treatment x PCE X experience effects

# Experienced loss

models_out_int_triple <- lapply(1:length(dv_vars_lag),function(i){
  tmp <- survey %>% 
    filter(term_status == "complete") %>%
    mutate(dv = scale(get(dv_vars_lag[i])),
           iv = as.numeric(scale(get(iv_vars_lag[i]))))
  
  tmp <- lm(dv ~ ideo_treat_binary * pce_treat_binary * experienced_loss + iv, data = tmp) %>% 
    margins::margins(variables = "ideo_treat_binary",
                     at = list(experienced_loss = c(T,F),
                               pce_treat_binary = c(T, F))) %>%
    summary() %>% tibble()  
  
  tmp$outcome <- dv_vars_lag[i]
  tmp
}
)
models_out_int_triple <- do.call("rbind", models_out_int_triple)

experienced_loss_plot <- models_out_int_triple %>%
  mutate(experienced_loss = ifelse(experienced_loss, "Experienced COVID Loss", "No COVID Loss"),
         experienced_loss = factor(experienced_loss, levels = c("No COVID Loss", "Experienced COVID Loss")),
         pce_treat_binary = ifelse(pce_treat_binary, "Primed", "Not primed"),
         pce_treat_binary = factor(pce_treat_binary, levels = c("Not primed", "Primed")),
         outcome = gsub("_post|_post_1","",outcome),
         outcome = factor(outcome, levels = c("redistSelf", "taxSpendSelf","jobsForAll","RoGIndex","redistIndex"))
         ) %>%
  ggplot(aes(x = AME, xmin = lower, xmax = upper, y = pce_treat_binary, col = experienced_loss)) + geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~outcome, ncol = 5) + 
  geom_vline(xintercept = 0, lty = 2) + 
  theme_bw() + ylab("") +
  scale_color_manual("", values = c("black", "gray")) + 
  xlab("Conditional treatment effect of ideological rhetoric") + 
  theme(legend.position = "bottom")


# Experienced support

models_out_int_triple <- lapply(1:length(dv_vars_lag),function(i){
  tmp <- survey %>% 
    filter(term_status == "complete") %>%
    mutate(dv = scale(get(dv_vars_lag[i])),
           iv = as.numeric(scale(get(iv_vars_lag[i]))))
  
  tmp <- lm(dv ~ ideo_treat_binary * pce_treat_binary * experienced_receipt + iv, data = tmp) %>% 
    margins::margins(variables = "ideo_treat_binary",
                     at = list(experienced_receipt = c(T,F),
                               pce_treat_binary = c(T, F))) %>%
    summary() %>% tibble()  
  
  tmp$outcome <- dv_vars_lag[i]
  tmp
}
)
models_out_int_triple <- do.call("rbind", models_out_int_triple)

experienced_receipt_plot <- models_out_int_triple %>%
  mutate(experienced_receipt = ifelse(experienced_receipt, "Experienced COVID Support", "No COVID Support"),
         experienced_receipt = factor(experienced_receipt, levels = c("No COVID Support", "Experienced COVID Support")),
         pce_treat_binary = ifelse(pce_treat_binary, "Primed", "Not primed"),
         pce_treat_binary = factor(pce_treat_binary, levels = c("Not primed", "Primed")),
         outcome = gsub("_post|_post_1","",outcome),
         outcome = factor(outcome, levels = c("redistSelf", "taxSpendSelf","jobsForAll","RoGIndex","redistIndex"))
  ) %>%
  ggplot(aes(x = AME, xmin = lower, xmax = upper, y = pce_treat_binary, col = experienced_receipt)) + geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~outcome, ncol = 5) + 
  geom_vline(xintercept = 0, lty = 2) + 
  theme_bw() + ylab("") +
  scale_color_manual("", values = c("black", "gray")) + 
  xlab("Conditional treatment effect of ideological rhetoric") + 
  theme(legend.position = "bottom")
  

pdf("../generated_images/treatment_triple_interaction_effects.pdf", width = 10, height = 6)
gridExtra::grid.arrange(experienced_loss_plot, experienced_receipt_plot, ncol = 1)
dev.off()

